/*

This file sets up tables

*/


**************** Make Sample 1 **************** 

use "/projects/data/dataSTATA/combined/combined_productivity_energy_4digit_4m.dta", clear

* drop obs that would be excluded from regressions --> generate constant sample
ivreghdfe lelec_int  (lelec_price lelec_price_init = instr_B_cl instr_B_ng instr_B_pa  instr_B_cl_init instr_B_ng_init instr_B_pa_init) [pw=wt_cmf], absorb(i.y##i.firstyear##i.bestnaics i.fipsst##i.bestnaics) cluster(fipsst)
keep if e(sample)

ivreghdfe lco2_int_i_both  (lelec_price lelec_price_init = instr_B_cl instr_B_ng instr_B_pa  instr_B_cl_init instr_B_ng_init instr_B_pa_init) [pw=wt_cmf], absorb(i.y##i.firstyear##i.bestnaics i.fipsst##i.bestnaics) cluster(fipsst)
keep if e(sample)

ivreghdfe lbtu_int_i_both  (lelec_price lelec_price_init = instr_B_cl instr_B_ng instr_B_pa  instr_B_cl_init instr_B_ng_init instr_B_pa_init) [pw=wt_cmf], absorb(i.y##i.firstyear##i.bestnaics i.fipsst##i.bestnaics) cluster(fipsst)
keep if e(sample)

reghdfe lelec_price_init instr_B_cl instr_B_ng instr_B_pa instr_B_cl_init instr_B_ng_init instr_B_pa_init   [pw=wt_cmf], absorb(i.y##i.firstyear##i.bestnaics i.fipsst##i.bestnaics) cluster(fipsst)
keep if e(sample)

gen age = year - firstyear

*>--------- define rounding function for regs---------------------<
* copied from rounding_regs.do Census file provided

capture program drop round
program define round, eclass
   mat A = e(b)
   forval j = 1/`=colsof(matrix(A))' { 
      mat A[1,`j'] = round(A[1,`j'],10^(floor(log10(abs(A[1,`j'])))-3))
   } 
   ereturn repost b = A

   if e(N)>=15&e(N)<=99 ereturn scalar N = round(e(N), 10)
   if e(N)>=100&e(N)<=999 ereturn scalar N = round(e(N), 50)
   if e(N)>=1000&e(N)<=9999 ereturn scalar N = round(e(N), 100)
   if e(N)>=10000&e(N)<=99999 ereturn scalar N = round(e(N), 500)
   if e(N)>=100000&e(N)<=999999 ereturn scalar N = round(e(N), 1000)
   if e(N)>=1000000 ereturn scalar N = round(e(N),10^(floor(log10(e(N)))-3))

   if e(N)<15 {
   scalar N_less15=1
   ereturn scalar N=.
   }

   if e(df_r)>=15&e(df_r)<=99 ereturn scalar df_r = round(e(df_r), 10)
   if e(df_r)>=100&e(df_r)<=999 ereturn scalar df_r = round(e(df_r), 50)
   if e(df_r)>=1000&e(df_r)<=9999 ereturn scalar df_r = round(e(df_r), 100)
   if e(df_r)>=10000&e(df_r)<=99999 ereturn scalar df_r = round(e(df_r), 500)
   if e(df_r)>=100000&e(df_r)<=999999 ereturn scalar df_r = round(e(df_r), 1000)
   if e(df_r)>=1000000 ereturn scalar df_r = round(e(df_r),10^(floor(log10(e(df_r)))-3))

   if e(df_r)<15 {
   scalar df_less15=1
   ereturn scalar df_r=.
   }

   capture ereturn scalar r2 = round(e(r2),10^(floor(log10(e(r2)))-3))
   capture ereturn scalar r2_a = round(e(r2_a),10^(floor(log10(e(r2_a)))-3))
   capture ereturn scalar F = round(e(F),10^(floor(log10(e(F)))-3))
   capture ereturn scalar rmse = round(e(rmse),10^(floor(log10(e(rmse)))-3))
   capture ereturn scalar mss = round(e(mss),10^(floor(log10(e(mss)))-3))
   capture ereturn scalar rss = round(e(rss),10^(floor(log10(e(rss)))-3))
   capture ereturn scalar rkf = round(e(rkf),10^(floor(log10(e(rkf)))-3))

end

**************** Add requested LIML regressions **************** 

********** Add panel to table 3 that runs LIML on panel B regs (IVregs) ******

ivreghdfe lelec_int (lelec_price lelec_price_init = instr_B_cl instr_B_ng instr_B_pa  instr_B_cl_init instr_B_ng_init instr_B_pa_init) [pw=wt_cmf], absorb(i.y##i.firstyear##i.bestnaics i.fipsst##i.bestnaics) cluster(fipsst) liml
eststo: round
ivreghdfe lelec_int (lelec_price lelec_price_init = instr_B_cl instr_B_ng instr_B_pa  instr_B_cl_init instr_B_ng_init instr_B_pa_init) [pw=wt_cmf], absorb(i.y##i.firstyear##i.bestnaics i.fipsst##i.bestnaics) cluster(fipsst) liml, if elec_perc_ind > 0.7
eststo: round
ivreghdfe lco2_int_i_both (lelec_price lelec_price_init = instr_B_cl instr_B_ng instr_B_pa  instr_B_cl_init instr_B_ng_init instr_B_pa_init) [pw=wt_cmf], absorb(i.y##i.firstyear##i.bestnaics i.fipsst##i.bestnaics) cluster(fipsst) liml
eststo: round
ivreghdfe lbtu_int_i_both (lelec_price lelec_price_init = instr_B_cl instr_B_ng instr_B_pa  instr_B_cl_init instr_B_ng_init instr_B_pa_init) [pw=wt_cmf], absorb(i.y##i.firstyear##i.bestnaics i.fipsst##i.bestnaics) cluster(fipsst) liml
eststo: round

esttab using "/projects/disclosure/20250205/output/robustness_energy_intensity_liml_T13T26.xls", b(3) se(3) obslast replace star(* 0.10 ** 0.05 *** 0.01) keep(lelec_price lelec_price_init) mtitles("Electricity Intensity" "Electricity Intensity, Elec. Intensive Industries" "CO2 Intensity" "BTU Intensity") title("LIML") coeflabels(lelec_price log(Current_Electricity_Price) lelec_price_init log(Initial_Electricity_Price)) scalars("rkf K-P F stat") sfmt(1)
eststo clear

********** Add panel to table 5 that runs LIML on panel B regs (IVregs) ********** 

ivreghdfe we_scale (lelec_price lelec_price_init = instr_B_cl instr_B_ng instr_B_pa  instr_B_cl_init instr_B_ng_init instr_B_pa_init) [pw=wt_cmf], absorb(i.y##i.firstyear##i.bestnaics i.fipsst##i.bestnaics) cluster(fipsst) liml
eststo: round
ivreghdfe we_scale (lelec_price lelec_price_init = instr_B_cl instr_B_ng instr_B_pa  instr_B_cl_init instr_B_ng_init instr_B_pa_init) [pw=wt_cmf], absorb(i.y##i.firstyear##i.bestnaics i.fipsst##i.bestnaics) cluster(fipsst) liml, if elec_perc_ind > 0.7
eststo: round
ivreghdfe wh (lelec_price lelec_price_init = instr_B_cl instr_B_ng instr_B_pa  instr_B_cl_init instr_B_ng_init instr_B_pa_init) [pw=wt_cmf], absorb(i.y##i.firstyear##i.bestnaics i.fipsst##i.bestnaics) cluster(fipsst) liml
eststo: round
ivreghdfe wh (lelec_price lelec_price_init = instr_B_cl instr_B_ng instr_B_pa  instr_B_cl_init instr_B_ng_init instr_B_pa_init) [pw=wt_cmf], absorb(i.y##i.firstyear##i.bestnaics i.fipsst##i.bestnaics) cluster(fipsst) liml, if elec_perc_ind > 0.7
eststo: round

esttab using "/projects/disclosure/20250205/output/robustness_productivity_liml_T13T26.xls", b(3) se(3) obslast replace star(* 0.10 ** 0.05 *** 0.01) keep(lelec_price lelec_price_init) mtitles("Energy Productivity" "Energy Productivity, Elec. Intensive Industries" "Total Factor Productivity" "Total Factor Productivity, Elec. Intensive Industries") title("LIML") coeflabels(lelec_price log(Current_Electricity_Price) lelec_price_init log(Initial_Electricity_Price)) scalars("rkf K-P F stat") sfmt(1)
eststo clear


**************** Add electricity shares from MECS sample **************** 


tempfile reg_sample
save `reg_sample', replace

local i = 0
foreach y of numlist 1985 1988 1991 1994 1998 2002 2006 2010 2014 {
	if inlist(`y', 2002) {
		import sas using "/data/economic/cmf/`y'/cmf`y'.sas7bdat", case(lower) clear
	}
	else {
		import sas using "/data/economic/asm/`y'/asm`y'.sas7bdat", case(lower) clear
	}
	replace year = `y'
	keep lbdnum lbdrevlnk year
	
	if `i' == 1 append using "/projects/data/dataSTATA/economic/asm_lbdnums.dta"
	local i = 1
	sa "/projects/data/dataSTATA/economic/asm_lbdnums.dta", replace
}

use "/projects/data/dataSTATA/economic/mecs_merged.dta", clear

tempfile mecs_firms
save `mecs_firms', replace

use "/projects/data/dataSTATA/economic/asm_lbdnums.dta", clear
drop if lbdnum==""

merge 1:m lbdrevlnk year using `mecs_firms'
keep if _merge==3
drop _merge
destring lbdnum, replace

merge 1:1 lbdnum year using `reg_sample'
keep if _merge==3


** Calculate share two ways
preserve
mat M = J(2, 2,.)


* electricity expenditures as fraction of tot energy expenditures ($)
do /projects/programs/codeSTATA/B_analysis/rounding_4sigdig_v2.do elec_btu btu_tot
gen elec_btu_sh1 = elec_btu / btu_tot
sum elec_btu_sh1
gen mean_elec_btu_sh1 = `r(mean)'
gen N_elec_btu_sh1 = `r(N)'
do /projects/programs/codeSTATA/B_analysis/rounding_4sigdig_v2.do mean_elec_btu_sh1 
do /projects/programs/codeSTATA/B_analysis/rounding_N.do N_elec_btu_sh1
sum mean_elec_btu_sh1 
mat M[1,1] = `r(mean)'
sum N_elec_btu_sh1
mat M[1, 2] = `r(mean)'


* electricity energy as fraction of tot energy consumption (BTU)

do /projects/programs/codeSTATA/B_analysis/rounding_4sigdig_v2.do pe btu_raw 

* Convert pe to BTUs - done in a seconary file for disclosure purposes
do projects/programs/codeSTATA/B_analysis/7_reducedform_disclosure_unitsconversion.do
gen elec_btu_sh2 = elec_btu_asm / (elec_btu_asm + btu_raw)
sum elec_btu_sh2

gen mean_elec_btu_sh2 = `r(mean)'
gen N_elec_btu_sh2 = `r(N)'
do /projects/programs/codeSTATA/B_analysis/rounding_4sigdig_v2 mean_elec_btu_sh2 
do /projects/programs/codeSTATA/B_analysis/rounding_N.do N_elec_btu_sh2 
sum mean_elec_btu_sh2
mat M[2,1] = `r(mean)'
sum N_elec_btu_sh2 
mat M[2,2] = `r(mean)'

matrix colnames M = Mean Obs
matrix rownames M = Elec_BTU_Share_MECs Elec_BTU_Share_ASM
esttab matrix(M) using "/projects/disclosure/20250205/output/energy_shares_T13T26.xls", replace

restore
